## Packages ####################################################################
library("coin")
library("plm")
library("lmtest")
library("car")
## library("memisc") loaded later
set.seed(0)


## Convenience functions #######################################################

## custom coefficient table
cftable <- function(obj, type = "single", var = "st") {
  str <- sprintf("%s%s:", var, type)
  cf <- coeftest(obj, vcov = vcovBK)
  ix <- grep(str, rownames(cf))
  if((type == "team") && ("twomaleyes" %in% rownames(cf))) ix <- c(ix, which(rownames(cf) == "twomaleyes"))
  rval <- cf[ix, , drop = FALSE]
  rownames(rval) <- gsub(str, "", rownames(rval), fixed = TRUE)
  class(rval) <- class(cf)
  attr(rval, "contrasts") <- obj$contrasts
  attr(rval, "xlevels") <- obj$xlevels
  attr(rval, "call") <- obj$call
  attr(rval, "nobs") <- as.numeric(rowSums(table(obj$model[[var]], attr(obj$model, "index")$teamid) > 0)[type])
  rval
}

## memisc::mtable plugin for custom coeftest table
getSummary.coeftest <- function(obj, alpha = 0.05, ...) {
  ## coefficient matrix and confidence interval
  cf <- cbind(obj, obj[,1] + qnorm(alpha/2) * obj[,2], obj[,1] + qnorm(1 - alpha/2) * obj[,2])
  colnames(cf) <- c("est", "se", "stat", "p", "lwr", "upr")

  ## further summary statistics
  sstat <- c("N" = attr(obj, "nobs"))

  ## return everything
  return(list(
    coef = cf,
    sumstat = sstat,
    contrasts = attr(obj, "contrasts"),
    xlevels = attr(obj, "xlevels"),
    call = attr(obj, "call")
  ))
}


## Data ########################################################################

## read CSV file
mla <- read.csv("MyopicLossAversion.csv")

## for convenience: abbreviate treatment variable names
names(mla) <- gsub("treatment_", "", names(mla), fixed = TRUE)

## add payouts in each round
for(i in 1:9) mla[[paste0("payout", i)]] <- (100 - mla[[paste0("invest", i)]]) +
  mla[[paste0("win", i)]] * 3.5 * mla[[paste0("invest", i)]]

## add transformed covariates
mla <- transform(mla,

  ## cumulative wins in all previous rounds
  cwin1 = 0,
  cwin2 = win1,
  cwin3 = win1 + win2,
  cwin4 = win1 + win2 + win3,
  cwin5 = win1 + win2 + win3 + win4,
  cwin6 = win1 + win2 + win3 + win4 + win5,
  cwin7 = win1 + win2 + win3 + win4 + win5 + win6,
  cwin8 = win1 + win2 + win3 + win4 + win5 + win6 + win7,
  cwin9 = win1 + win2 + win3 + win4 + win5 + win6 + win7 + win8,

  ## lagged wins in last previous round
  lwin1 = 0,
  lwin2 = win1,
  lwin3 = win2,
  lwin4 = win3,
  lwin5 = win4,
  lwin6 = win5,
  lwin7 = win6,
  lwin8 = win7,
  lwin9 = win8,
  
  ## wins in last block
  bwin1 = 0,
  bwin2 = 0,
  bwin3 = 0,
  bwin4 = win1 + win2 + win3,
  bwin5 = 0,
  bwin6 = 0,
  bwin7 = win4 + win5 + win6,
  bwin8 = 0,
  bwin9 = 0,

  ## cumulative wins in all previous rounds
  cpayout1 = 0,
  cpayout2 = payout1,
  cpayout3 = payout1 + payout2,
  cpayout4 = payout1 + payout2 + payout3,
  cpayout5 = payout1 + payout2 + payout3 + payout4,
  cpayout6 = payout1 + payout2 + payout3 + payout4 + payout5,
  cpayout7 = payout1 + payout2 + payout3 + payout4 + payout5 + payout6,
  cpayout8 = payout1 + payout2 + payout3 + payout4 + payout5 + payout6 + payout7,
  cpayout9 = payout1 + payout2 + payout3 + payout4 + payout5 + payout6 + payout7 + payout8,

  ## lagged payouts in last previous round
  lpayout1 = 0,
  lpayout2 = payout1,
  lpayout3 = payout2,
  lpayout4 = payout3,
  lpayout5 = payout4,
  lpayout6 = payout5,
  lpayout7 = payout6,
  lpayout8 = payout7,
  lpayout9 = payout8,
  
  ## payouts in last block
  bpayout1 = 0,
  bpayout2 = 0,
  bpayout3 = 0,
  bpayout4 = payout1 + payout2 + payout3,
  bpayout5 = 0,
  bpayout6 = 0,
  bpayout7 = payout4 + payout5 + payout6,
  bpayout8 = 0,
  bpayout9 = 0,

  ## (re-)level categorical covariates
  gender   = factor(gender, levels = c("male", "female")),
  school   = factor(school),
  yo       = factor(yo, levels = c("Younger", "Older")),

  ## age centered: approximately within age group
  age_yo   = age - c(17, 12)[yo]
)

## additional IDs for single players (to facilitate subsequent aggregation)
mla$teamid[is.na(mla$teamid)] <- 1:sum(is.na(mla$teamid)) + max(mla$teamid, na.rm = TRUE)
mla$teamid <- factor(mla$teamid)

## aggregate within teams
## categorical variables
amla <- with(mla, data.frame(
  teamid = factor(tapply(teamid, teamid, unique), labels = levels(teamid)),
  lh = factor(tapply(lh, teamid, unique), labels = levels(lh)),
  yo = factor(tapply(yo, teamid, unique), labels = levels(yo)),
  st = factor(tapply(st, teamid, unique), labels = levels(st)),
  school = factor(tapply(school, teamid, unique), labels = levels(school)),
  gender = factor(tapply(gender, teamid, function(x) if(length(x) > 1 && x[1] != x[2]) "female/male" else as.character(x)[1])),
  male = factor(tapply(gender, teamid, function(x) "male" %in% x), levels = c(FALSE, TRUE), labels = c("no", "yes")),
  twomale = factor(tapply(gender, teamid, function(x) if(length(x) > 1 && x[1] == x[2] && x[1] == "male") "yes" else "no"))
))
## numeric variables
for(n in names(mla)[which(sapply(mla, is.numeric))]) amla[[n]] <- as.vector(tapply(mla[[n]], mla[["teamid"]], mean))


## reshape into "long" format
var <- grep("[1-9]", names(amla))
vnm <- names(amla)[var]
vnm <- gsub("[1-9]", "", vnm)
lmla <- reshape(amla, direction = "long",
  idvar = "teamid", timevar = "round",
  varying = split(var, factor(vnm, levels = unique(vnm))),
  v.names = unique(vnm))
lmla <- lmla[order(lmla$teamid),]
lmla$age <- lmla$age_yo


## for subjects in LONG treatment use lwin = bwin / lpayout = bpayout
lmla$lwin[lmla$lh == "LONG"] <- lmla$bwin[lmla$lh == "LONG"] / 3
lmla$lpayout[lmla$lh == "LONG"] <- lmla$bpayout[lmla$lh == "LONG"] / 3
## omit LONG investment decision in non-decision rounds
lmla <- subset(lmla, lh == "SHORT" | round %in% c(1, 4, 7))
lmla$round <- as.numeric(lmla$round)


## panel data _with_ first round
pmla <- pdata.frame(lmla, c("teamid", "round"))


## and omitting first round
pmla1 <- pdata.frame(subset(lmla, round > 1), c("teamid", "round"))


## table needed later for Figure 2 (but needs to be computed before attaching memisc)
nd2e <- aggregate(invest ~ st + yo + male, data = pmla, FUN = mean)


## Table 1 #####################################################################
ftable(xtabs(~ yo + st + lh, data = mla)[, ,2:1])


## Table 2 #####################################################################
tab2 <- with(amla, tapply(minvest, interaction(lh, st, yo), function(x)
  c("Mean" = mean(x), "Median" = median(x), "St. dev." = sqrt(var(x)/length(x)), "# units" = length(x))))
tab2 <- do.call("cbind", tab2)
round(tab2, digits = 2)


## Nonparametric inference (Section 3.1) #######################################

## first analysis of myopic loss aversion using nonparametric tests
wilcox_test(minvest ~ lh, data = amla, subset = amla$yo == "Younger" & amla$st == "single")
wilcox_test(minvest ~ lh, data = amla, subset = amla$yo == "Older"   & amla$st == "single")
wilcox_test(minvest ~ lh, data = amla, subset = amla$yo == "Younger" & amla$st == "team")
wilcox_test(minvest ~ lh, data = amla, subset = amla$yo == "Older"   & amla$st == "team")

## and analysis of team effect
wilcox_test(minvest ~ st, data = amla, subset = amla$yo == "Younger" & amla$lh == "SHORT")
wilcox_test(minvest ~ st, data = amla, subset = amla$yo == "Older"   & amla$lh == "SHORT")
wilcox_test(minvest ~ st, data = amla, subset = amla$yo == "Younger" & amla$lh == "LONG")
wilcox_test(minvest ~ st, data = amla, subset = amla$yo == "Older"   & amla$lh == "LONG")


## Figure 1 ######################################

## empirical quantities from Table 2
mn <- with(amla, tapply(minvest, list(st, lh, yo), mean))[,2:1,]
md <- with(amla, tapply(minvest, list(st, lh, yo), median))[,2:1,]
se <- with(amla, tapply(minvest, list(st, lh, yo), function(x) sqrt(var(x)/length(x))))[,2:1,]

## payouts (me, se): overall, by yo, by st * yo * lh
meanse <- function(x, digits = 1) paste(format(round(mean(x), digits = digits), nsmall = digits),
  " (", format(round(sd(x), digits = digits), nsmall = digits), ")", sep = "")
amla$payout <- rowSums(amla[, paste("payout", 1:9, sep = "")])
pay1 <- with(amla, meanse(payout))
pay2 <- with(amla, tapply(payout, yo, meanse))
pay8 <- with(amla, tapply(payout, list(st, lh, yo), meanse))[,2:1,]
ftable(pay8, col.vars = 1)

# rb <- hcl(c(0, 260), 40, 70, fixup = FALSE)
rb <- gray(c(0.75, 0.35))
par(mfrow = c(1, 2), mar = c(3, 0, 3, 0), oma = c(0, 4, 0, 3), yaxs = "i")
plot(0, 0, xlim = c(0.8, 6.2), ylim = c(0, 80), type = "n", axes = FALSE, xlab = "", ylab = "")
abline(h = 1:8 * 10, col = "lightgray")
xat <- barplot(mn[,,1], beside = TRUE, ylim = c(0, 80),
  main = "Younger group", names = c("SHORT", "LONG"),
  col = rb, add = TRUE, axes = FALSE)
for(i in 1:2) { for(j in 1:2) {
  arrows(xat[i,j], mn[i,j,1] - se[i,j,1],
    xat[i,j], mn[i,j,1] + se[i,j,1],
    angle = 90, code = 3, length = 0.05)
}}
axis(2)
box()
mtext("Average amount invested", side = 2, line = 3)
legend("topleft", c("Individual", "Team"), pch = 22, pt.bg = rb, pt.cex = 2.5, bty = "n")
plot(0, 0, xlim = c(0.8, 6.2), ylim = c(0, 80), type = "n", axes = FALSE, xlab = "", ylab = "")
abline(h = 1:8 * 10, col = "lightgray")
xat <- barplot(mn[,,2], beside = TRUE, ylim = c(0, 80),
  main = "Older group", names = c("SHORT", "LONG"),
  col = rb, add = TRUE, axes = FALSE)
for(i in 1:2) { for(j in 1:2) {
  arrows(xat[i,j], mn[i,j,2] - se[i,j,2],
    xat[i,j], mn[i,j,2] + se[i,j,2],
    angle = 90, code = 3, length = 0.05)
}}
axis(4)
box()


## Table 3 #####################################################################

## separate regressions by st (== interactions of st with all regressors except twomale dummy and school fixed effects)
tab3_m1 <- plm(invest ~ st / (yo * lh + male + yo * age) + twomale + school, data = pmla, model = "pooling")
tab3_m12<- plm(invest ~ st / (male + yo * age)           + twomale + school, data = pmla, model = "pooling")
tab3_m2 <- plm(invest ~ st / (male + yo * age)                     + school, data = pmla, model = "pooling")
waldtest(tab3_m1, tab3_m12, tab3_m2, vcov = vcovBK)

## model table (via memisc)
library("memisc")
setSummaryTemplate("coeftest" = c(
  "N" = "($N:d)"
))

## Table 3 (note: memisc uses other significance star coding!)
mtable(
  "1_single" = cftable(tab3_m1, "single"),
  "1_team" = cftable(tab3_m1, "team"),
  "2_single" = cftable(tab3_m2, "single"),
  "2_team" = cftable(tab3_m2, "team")
)


## Figure 2 ####################################################################
## model-based investment
tab3_m2_lm <- lm(terms(tab3_m2), data = pmla)
nd2 <- expand.grid(
  male = levels(pmla$male),
  st = levels(pmla$st),
  yo = levels(pmla$yo),
  school = "4"
)
nd2$age <- with(pmla, tapply(age, yo, mean))[nd2$yo]
nd2$invest <- predict(tab3_m2_lm, nd2)
## exploratory version (already computed above)
# nd2e <- aggregate(invest ~ st + yo + male, data = pmla, FUN = mean)

par(mfrow = c(1, 2), mar = c(3, 0, 3, 0), oma = c(0, 4, 0, 3), yaxs = "r")
rb <- hcl(c(260, 260, 0, 0), c(0, 0, 0, 0), c(15, 60, 25, 70), fixup = FALSE)
gr <- hcl(0, 0, 100)
plot(invest ~ as.numeric(st), data = nd2e, type = "n",
  main = "Younger group", axes = FALSE, xlim = c(0.5, 2.5), ylim = c(40, 75),
  xlab = "", ylab = "")
abline(h = seq(40, 75, by = 5), col = "lightgray")
lines(invest ~ as.numeric(st), subset = yo == "Younger" & male == "yes",
  data = nd2e, type = "b", pch = 17, col = rb[2], lty = 2, lwd = 1.5)
lines(invest ~ as.numeric(st), subset = yo == "Younger" & male == "no",
  data = nd2e, type = "b", pch = 19, col = rb[4], lty = 2, lwd = 1.5)
lines(invest ~ as.numeric(st), subset = yo == "Younger" & male == "yes",
  data = nd2, type = "b", pch = 17, col = rb[1], lty = 1, lwd = 1.5)
lines(invest ~ as.numeric(st), subset = yo == "Younger" & male == "no",
  data = nd2, type = "b", pch = 19, col = rb[3], lty = 1, lwd = 1.5)
axis(1, at = 1:2, labels = c("Individual", "Team"))
axis(2)
box()
mtext("Amount invested", side = 2, line = 3)
legend("topleft", c("Male, predicted", "Male, observed", "Female, predicted", "Female, observed"),
  col = rb, lty = c(1, 2, 1, 2), pch = c(17, 17, 19, 19), lwd = 1.5, bty = "n")
plot(invest ~ as.numeric(st), data = nd2e, type = "n", main = "Older group", axes = FALSE,
  xlim = c(0.5, 2.5), ylim = c(40, 75), xlab = "", ylab = "")
abline(h = seq(40, 75, by = 5), col = "lightgray")
lines(invest ~ as.numeric(st), subset = yo == "Older" & male == "yes",
  data = nd2e, type = "b", pch = 17, col = rb[2], lty = 2, lwd = 1.5)
lines(invest ~ as.numeric(st), subset = yo == "Older" & male == "no",
  data = nd2e, type = "b", pch = 19, col = rb[4], lty = 2, lwd = 1.5)
lines(invest ~ as.numeric(st), subset = yo == "Older" & male == "yes",
  data = nd2, type = "b", pch = 17, col = rb[1], lty = 1, lwd = 1.5)
lines(invest ~ as.numeric(st), subset = yo == "Older" & male == "no",
  data = nd2, type = "b", pch = 19, col = rb[3], lty = 1, lwd = 1.5)
axis(1, at = 1:2, labels = c("Individual", "Team"))
axis(4)
box()


## Table 4 #####################################################################

tab4_m1 <- plm(invest ~ st / (male + yo * age)                + school, data = pmla1, model = "pooling")
tab4_m2 <- plm(invest ~ st / (male + yo * age + lwin + cwin ) + school, data = pmla1, model = "pooling")
waldtest(tab4_m1, tab4_m2, vcov = vcovBK)

mtable(
  "1_single" = cftable(tab4_m1, "single"),
  "1_team" = cftable(tab4_m1, "team"),
  "2_single" = cftable(tab4_m2, "single"),
  "2_team" = cftable(tab4_m2, "team")
)

# test for differences in short- and long-sighted effects
tab4_m2_lh <- plm(invest ~ st / (male + yo * age + lh * (lwin + cwin))   + school, data = pmla1, model = "pooling")
waldtest(tab4_m2, tab4_m2_lh, vcov = vcovBK)
tab4_m2_ge <- plm(invest ~ st / (male + yo * age + male * (lwin + cwin)) + school, data = pmla1, model = "pooling")
waldtest(tab4_m2, tab4_m2_ge, vcov = vcovBK)
tab4_m2_yo <- plm(invest ~ st / (male + yo * age + yo * (lwin + cwin))   + school, data = pmla1, model = "pooling")
waldtest(tab4_m2, tab4_m2_yo, vcov = vcovBK)


## Figure 3 ####################################################################

tab4_m3_lm <- lm(terms(tab4_m2), data = pmla1)
nd3 <- expand.grid(
  lwin = 0:1,
  cwin = c(1, 3),
  st = levels(pmla1$st),
  yo = levels(pmla1$yo),
  school = "4"
)
nd3$male <- factor(rep("yes", nrow(nd3)), levels = c("no", "yes"))
nd3$age <- with(pmla1, tapply(age, yo, mean))[nd3$yo]
nd3$invest <- predict(tab4_m3_lm, nd3)

par(mfrow = c(1, 2), mar = c(3, 0, 3, 0), oma = c(0, 4, 0, 3))
# rb <- hcl(c(0, 260), 90, 40, fixup = FALSE)
rb <- hcl(c(0, 260), 0, c(15, 45), fixup = FALSE)
plot(invest ~ as.numeric(st), data = nd3, type = "n",
  main = "Younger group", axes = FALSE, xlim = c(0.5, 2.5), ylim = c(44, 82),
  xlab = "", ylab = "")
abline(h = seq(40, 85, by = 5), col = "lightgray")
lines(invest ~ as.numeric(st), subset = yo == "Younger" & lwin == 0 & cwin <= 2,
  data = nd3, type = "b", pch = 17, col = rb[1], lty = 2, lwd = 1.5)
lines(invest ~ as.numeric(st), subset = yo == "Younger" & lwin == 0 & cwin >  2,
  data = nd3, type = "b", pch = 17, col = rb[1], lty = 1, lwd = 1.5)
lines(invest ~ as.numeric(st), subset = yo == "Younger" & lwin == 1 & cwin <= 2,
  data = nd3, type = "b", pch = 19, col = rb[2], lty = 2, lwd = 1.5)
lines(invest ~ as.numeric(st), subset = yo == "Younger" & lwin == 1 & cwin >  2,
  data = nd3, type = "b", pch = 19, col = rb[2], lty = 1, lwd = 1.5)
axis(1, at = 1:2, labels = c("Individual", "Team"))
axis(2)
box()
mtext("Amount invested", side = 2, line = 3)
legend("topleft", c("Short-sight = 0, Long-sight = 3", "Short-sight = 0, Long-sight = 1", "Short-sight = 1, Long-sight = 3", "Short-sight = 1, Long-sight = 1"),
  col = rb[c(1, 1, 2, 2)], lty = c(1, 2, 1, 2), pch = c(17, 17, 19, 19), lwd = 1.5, bty = "n")
plot(invest ~ as.numeric(st), data = nd3, type = "n",
  main = "Older group", axes = FALSE, xlim = c(0.5, 2.5), ylim = c(44, 82),
  xlab = "", ylab = "")
abline(h = seq(40, 85, by = 5), col = "lightgray")
lines(invest ~ as.numeric(st), subset = yo == "Older" & lwin == 0 & cwin <= 2,
  data = nd3, type = "b", pch = 17, col = rb[1], lty = 2, lwd = 1.5)
lines(invest ~ as.numeric(st), subset = yo == "Older" & lwin == 0 & cwin >  2,
  data = nd3, type = "b", pch = 17, col = rb[1], lty = 1, lwd = 1.5)
lines(invest ~ as.numeric(st), subset = yo == "Older" & lwin == 1 & cwin <= 2,
  data = nd3, type = "b", pch = 19, col = rb[2], lty = 2, lwd = 1.5)
lines(invest ~ as.numeric(st), subset = yo == "Older" & lwin == 1 & cwin >  2,
  data = nd3, type = "b", pch = 19, col = rb[2], lty = 1, lwd = 1.5)
axis(1, at = 1:2, labels = c("Individual", "Team"))
axis(4)
box()


## Footnote 3 ##################################################################

## subset with pocket money > 0 and add centered age (within grade)
mla0 <- subset(mla, pocketmoney > 0)
mla0 <- transform(mla0,
  age_gr   = age - tapply(age, grade, mean)[factor(grade)]
)

## pocket money model
pm <- lm(log(pocketmoney) ~ grade + age_gr, data = mla0)

## prediction different grades at mean age
pm_pred <- data.frame(grade = c(8, 10), age_gr = 0)
pm_pred$pocketmoney <- exp(predict(pm, newdata = pm_pred))
pm_pred
